require(mvtnorm)
require(OpenMx)
require(psych)
require(MASS)

inter <- 0                                                  # Generate the interaction parameter. We used 0, .5, 1, 1.5 & 2
nsim <- 10000                                               # Set the number of simulations

AllMZdata <- list(NULL)
AllDZdata <- list(NULL)

MZdatasets <- list(NULL)
DZdatasets <- list(NULL)


nv <- 1                                                     # Set the number of variables that you will be using (per person)
Vars <- c('t')                                              # Name the variable
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="")        # Create an object that will select the variables
ntv <- nv*2                                                 # Set the total number of variable that you will use (1 var * 2 twins)


mzmat <- matrix(c(1,1,0,0,0,0,
                  1,1,0,0,0,0,
                  0,0,1,1,0,0,
                  0,0,1,1,0,0,
                  0,0,0,0,1,0,
                  0,0,0,0,0,1),6)                           # Create a matrix of correlations for the genetic and environmental variance components for each twin (MZ)

dzmat <- matrix(c(1,1/2,0,0,0,0,
                  1/2,1,0,0,0,0,
                  0,0,1,1,0,0,
                  0,0,1,1,0,0,
                  0,0,0,0,1,0,
                  0,0,0,0,0,1),6)                           # Create a matrix of correlations for the genetic and environmental variance components for each twin (DZ)

Avar <- c(0:8/10,0:7/10,0:6/10,0:5/10)                      # Create an object with each of the possible Additive Genetic Parameters
Cvar <- c(8:0/10,7:0/10,6:0/10,5:0/10)                      # Create an object with each of the possible Common Environmental Parameters
Evar <- c(rep(.2,9),rep(.3,8),rep(.4,7),rep(.5,6))          # Create an object with each of the possible Unique Environmental Parameters


for (param in 1:30){
	for(sim in 1:nsim){
		mzSim  <- mvrnorm(n=1000, c(0,0,0,0,0,0), mzmat , empirical=T)     # Simulate genetic and environmental factor scores for each MZ twin
		dzSim  <- mvrnorm(n=1000, c(0,0,0,0,0,0), dzmat , empirical=T)     # Simulate genetic and environmental factor scores for each DZ twin

		VarCompA <- Avar[param]                                            # Set the Genetic Variance Component
		VarCompC <- Cvar[param]                                            # Set the Common Environmental Variance Component
		VarCompE <- Evar[param]                                            # Set the Unique Environmental Variance Component

		aPath <- sqrt(VarCompA)                                            # Set the Genetic Variance Path Coefficeint
		cPath <- sqrt(VarCompC)                                            # Set the Common Environmental Path Coefficeint
		ePath <- sqrt(VarCompE)                                            # Set the Unique Environmental Path Coefficeint

		MZdata <- cbind((mzSim[,1]*aPath+
		                 mzSim[,3]*cPath+
		                 mzSim[,5]*ePath+
		                 mzSim[,1]*mzSim[,5]*inter),
		                (mzSim[,2]*aPath+
		                 mzSim[,4]*cPath+
		                 mzSim[,6]*ePath+
		                 mzSim[,2]*mzSim[,6]*inter))                      # Construct the phenotype by adding the genetic/environmental factor score * the path + the bias due to the interaction for MZ twins

		DZdata <- cbind((dzSim[,1]*aPath+
		                 dzSim[,3]*cPath+
		                 dzSim[,5]*ePath+
		                 dzSim[,1]*dzSim[,5]*inter),
		                (dzSim[,2]*aPath+
		                 dzSim[,4]*cPath+
		                 dzSim[,6]*ePath+
		                 dzSim[,2]*dzSim[,6]*inter))                      # Construct the phenotype by adding the genetic/environmental factor score * the path + the bias due to the interaction for DZ twins

		colnames (MZdata) <- selVars                                      # Give the datasets column names
		colnames (DZdata) <- selVars
		
		
		AllMZdata[[sim]] <- MZdata	
		AllDZdata[[sim]] <- DZdata	
		
	}
	MZdatasets[[param]]<- AllMZdata
	DZdatasets[[param]]<- AllDZdata
}

#MZdatasets
#DZdatasets